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We introduce a new continuous cellular automaton that 
presents self-organized criticality. It is one-dimensional, to- 
tally deterministic, without any embedded randomness, not 
even in the initial conditions. This system is in the same uni- 
versality class as the Oslo rice pile, boundary driven interface 
depinning and the train model for earthquakes. Although the 
system is chaotic, in the thermodynamic limit chaos occurs 
only in a microscopic level. 



In 1987, Bak, Tang and Wiesenfeld showed that frac- 
tal behavior, that is, power-law distributions, can be ob- 
served in simple dissipative systems with many degrees 
of freedom without fine tuning of parameters JjJ . They 
called this phenomenon self-organized criticality (SOC). 
Until then, the studies of fractal structures were basically 
related to equilibrium systems where fractality appears 
only at special parameter values where a phase transition 
takes place. 

Since the pioneering work of Bak et. al, an enormous 
amount of numerical, theoretical and experimental stud- 
ies have been done in systems that present SOC. One of 
the most interesting experimental studies demonstrating 
the existence of SOC in Nature was done in a quasi-one- 
dimensional pile of rice by Frette et. al Q . They found 
that the occurrence of SOC depends on the shape of the 
rice. Only with sufficient elongated grains, avalanches 
with a power-law distribution occurred. If the rice had 
little asymmetry, a distribution described by a stretched 
exponential was seen. Christensen et. al ]3| introduced a 
model for the rice pile experiment in which the local crit- 
ical slope varies randomly between 1 and 2. They found 
that their model, known as the Oslo rice pile model, re- 
produced well the experimental results. 

A good understanding of the Oslo system was achieved 
by Paczuski and Boettcher B . They showed that it could 
be mapped exactly to a model for interface depinning 
where the interface is slowly pulled at one end through 
a medium with quenched random pinning forces. They 
found that the height of the interface maps to the num- 
ber of toppling events in the rice pile model. The crit- 
ical exponents of the two models were identical (within 
the error bars), showing that they were in the same uni- 
versality class. Paczuski and Boettcher also conjectured 
that the train model for earthquakes, which was intro- 
duced by Burridge and Knopoff @] , and studied in detail 
in is also in that same universality class. The train 
model is the only model that we know (besides the one 



we introduce here) that presents SOC and has no kind of 
embedded randomness. However, it is governed by cou- 
pled ordinary equations (ODE's), what makes its study 
very time consuming. 

A way of making a system governed by ODE's more 
amenable to computer simulations is to discretize it in 
time. This was done by Olami, Feder and Christensen 
(OFC) who introduced a continuous cellular automa- 
ton (CCA) to study the two-dimensional version of an- 
other Burridge and Knopoff model for earthquakes || . [A 
continuous cellular automaton in SOC is known in chaos 
theory as coupled lattice maps. These systems are char- 
acterized as having space and time variables defined in 
the domain of real and integer numbers, respectively.] In 
the OFC model, SOC is seen only in systems that have 
a geometry with dimensionality of at least two. That 
system is a variation of a model introduced by Nakanishi 
which has a one-dimensional geometry. However, the 
model introduced by Nakanishi does not present SOC, 
since the power-law distribution it presents has an upper 
cutoff that is unrelated to the system size. 

Here we introduce a new self-organized critical system, 
that is governed by a CCA (that is, the space is contin- 
uous and the time is discrete). It is one-dimensional and 
has no embedded randomness, not even in the initial con- 
ditions. We will show that our model belongs to the same 
universality class as the Oslo rice pile, boundary driven 
interface depinning and the train model. The importance 
of our results comes from the fact that we show that it is 
possible to map stochastic SOC systems to simple, dis- 
crete, chaotic systems, in which no randomness exists. 
Such an equivalence of a deterministic model with no 
embedded randomness which is chaotic with a stochastic 
model also occurs between the deterministic Kuramoto- 
Shivashinsky [lC| ] equation and the Langevin equation 
proposed by Kardar, Parisi and Zhang pd| ]. In our opin- 
ion, the train model governed by ODE's already achieved 
this H . However, because its simulation is very time con- 
suming, it will probably be impossible to find such equiv- 
alences for higher dimensional (two, three, etc.) systems, 
unless the system is discretized in time, as we do here. In 
fact, we are unaware of any studies on train-like systems 
with dimensionality higher than one. To the best of our 
knowledge, our model is the only SOC system introduced 
so far that is one-dimensional, totally deterministic and 
with discrete time. 

Another important result of this paper concerns to the 
fact that we show that although chaos exists in the model, 
it decreases as the system size increases, and in thermo- 
dynamic limit it exists only in a microscopic level. Con- 
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sequently, our results indicate that the fractal structures 
seen in Nature and supposed associated with SOC, could 
in principle result only from nonlinearities in those sys- 
tems, without any need for the presence of random irreg- 
ularities. Such fundamental questions are also found in 
equilibrium statistical mechanics, where it is unknown if 
probability theory is only a tool to describe phenomena 
that in principle could have been originated solely from 
microscopic chaos p|. 

The train model is shown in Fig. [l](a). It consists of a 
chain of blocks connect by harmonic springs. The blocks 
are on a rough surface with friction, and the first block is 
pulled slowly with a constant velocity by a driver. The 
dynamics of the model is as follows: suppose that all the 
blocks are initially at rest. As the driver pulls the first 
block, it remains stuck until the elastic force overcomes 
the static friction. When this occurs, the first block will 
move a little and will be stopped again by friction. Such 
small events (or earthquakes) will continue, and will in- 
crease the elastic force on the second block. There will be 
a moment when the elastic force on the second block will 
overcome the friction force, and then we will see an event 
involving two blocks. This kind of dynamics will continue 
with events involving three, four, or all the blocks in the 
system. 

The elastic force in block i is given by fi = Xi—% — 
2xi + Xi+\, where xi is the displacement of block i with 
respect to its equilibrium position (without losing gener- 
ality, the spring constant can be taken as equal to 1 H). 
The boundary conditions are xq = and xl+i = %l- 
After an earthquake, in which block i was displaced by 
Ax, the elastic forces in block i and in its nearest neigh- 
bors will be f! = /< - A/, and f' i±1 = f i±1 + A//2, 
respectively, where A/ = 2Ax. Thus, the force that 
is relaxed in block i is redistributed equally to its near 
nearest neighbors, implying conservation of elastic forces. 
This is embedded into the geometry of the system. How- 
ever, the model does dissipate energy through friction 
between the rough surface and the blocks. Consequently, 
the model is conservative with respect to elastic forces, 
and dissipative with respect to energy. This is one of 
the main distinctions between the train model and the 
other Burridge and Knopoff earthquake model studied 
by Nakanishi and OFC, in which neither the energy nor 
the forces are conserved. 

In the discretized version of the train model that we 
introduce here, each block i is characterized by a vari- 
able fi, which we will call force, with i = 1, L, and L 
being the number of blocks in the system. The boundary 
conditions are the same as the ones in the train model, 
which are given above. The dynamical evolution of the 
system is determined by the following algorithm: 

(1) Start the system by defining initial values for the 
variables which can be the same for all the blocks, so 
the they are below a chosen, fixed, threshold f t h- 

(2) Update the force in the first block by incrementing 



its force to the threshold value plus a fixed small value 
Sf, i.e., /i = f th + Sf (an event is going to be triggered). 

(3) Check the forces in each block. If a block i has 
fi > fth, update fi according to f[ = cj>(fi - f th ), where 
<f> is a given nonlinear function that has a parameter a. 
Increase the forces in its two nearest neighboring blocks 
according to f' l±1 = f l±1 + A//2, where A/ = fi - f[. 

(4) If f[ < f th for all the blocks, go to step (2) (the 
event has finished). Otherwise, go to step (3) (the event 
is still evolving). 

One can use either parallel or sequential update in the 
evolution of the system. We have verified that the critical 
exponents of the model do not depend on the type of 
update chosen. The systems is governed by L variables 
and has two parameters, a and Sf, since without losing 
generality we can take fth = 1 ■ The force in our model is 
supposed to mimic the combination of two forces in the 
train model, that is, the elastic and the friction forces. 
The elastic force is periodic, whereas the friction force 
in simulations is generally assumed to decrease with the 
velocity of the block. We have found numerically that 
<f> mimics the combination of these two force when it is 
a periodic function, since only in this way the system 
presents SOC behavior. So, the periodicity of the elastic 
force dominates over the form of the friction force. The 
models introduced by Nakanishi and OFC assume that 
is a strictly nonincreasing function. We have found 
that if we use a strictly nonincreasing function for <j), 
such as the one used in we observe in our model 
the presence of stretched exponentials, instead of power- 
laws Q. It is worth noting that in this situation our 
model reproduces the results of the distributions found 
with nearly round rice in which were also governed 
by stretched exponentials. 

The functional form we chose for <f> is shown in 
Fig. 0(b), which is given by, 4>(x) = 1 — a[x] where [x] 
denotes x modulo 1/a, that is, a sawtooth function. How- 
ever, we have tested several other periodic functions, and 
found that the SOC behavior we show here remains, that 
is, the results are robust, the essential ingredient being 
periodicity (not necessary a perfect one) for <j). 

In Fig. ^| we show the distribution of events involving s 
update steps, which is the size of the event, using parallel 
update. The events that involve all the blocks of the 
chain have been excluded from our analysis, since they do 
not belong to the same distribution, as expected. Before 
we start to compute the statistics, we wait until the last 
block has moved, in order to neglect transient effects. 
In (a) we show three different cases, with L = 512, in 
which we have varied one parameter at a time. The solid 
line refers to a = 4 and Sf = 0.1, the dashed line is for 
a = 4 and Sf = 0.01, and the short-dashed line refers 
to a = 2.5 and Sf = 0.01. We see that the small events 
have their own distribution, like in the Oslo model for 
rice piles ||. A careful analysis has shown that these 
small events have an exponential distribution M. As s 
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is increased the distribution becomes a power-law, which 
has a cutoff related to finite size effects, only. We find 
that the slope of the power law is independent of a and 
Sf. However, the crossover point s* of the exponential 
behavior to the power law one depends on a, but not 
on Sf. The frequency of the small events is inversely 
proportional to Sf. 

In Fig. ^(b) we show simulations keeping a and Sf 
fixed, with a = 4 and Sf — 0.1, and varying L. We see 
that increasing L, the range of the power-law increases. 
To fit the data to a single curve, we try the finite size scal- 
ing ansatz P(s, L) = s~ T G(s/L D ), where D and r are the 
basic exponents of the model |l3f| , defining its universality 
class. D and r are called the dimension and the size dis- 
tribution exponents, respectively. In our model we find 
that < s >~ L which results inr = 2— 1/D, The best 
fit for P(s,L) is found for t w 1.54 and D w 2.20. The 
data collapse for these values of the exponents is shown in 
Fig. ||(c). Within the error bars, these exponents are the 
same ones of the Oslo rice pile |||| , driven boundary in- 
terface depinning J| and the train model for earthquake 
H,|). Consequently, all these models, including the one 
we introduce here, are in the same universality class. 

In Fig. ||(a) we show the frequency of the events P(T) 
as a function of the its time duration T, for differ- 
ent system sizes. The parameters are the same as in 
Fig. ||(b). Again, we see a power-law distribution, except 
for the smallest events. A data collapse for the func- 
tion P(T, L) = T-yf T {T/L a ) is shown in Fig. |(b), with 
y = 1.84 and a = 1.40. These are the same exponents 
found in the Oslo rice pile. From conservation of proba- 
bility it follows that a(y — 1) = D(t — 1) in good agree- 
ment with our results. The results shown in Fig. 3 are, 
again, for parallel update, since in sequential update the 
time duration and event size coincide resulting in y = t 
and a = D. 

Using the method introduced by Benettin et al. [flit , 
we have calculated the largest Liapunov exponent, Ai, 
and the second largest Liapunov exponent, A2, of the 
system. If Ai is greater than zero, it implies that the 
system has a strong sensitivity to the initial conditions, 
and by definition, it is called chaotic. To study the Lia- 
punov exponents we have chosen sequential update. The 
reason for this is that the calculation of the Liapunov ex- 
ponent assumes, by its own definition, that small changes 
happen in the system in the time unit, and this is more 
consistent with sequential rather than with parallel up- 
date. In Fig. ||(a) and Fig. |](b) we show Ai and A 2 as 
a function of a for Sf = 0.1, and as a function of Sf for 
a = 4, respectively. In both cases L — 64. We have found 
that for a < 1 the system is in continuous motion, and 
therefore being impossible to define earthquakes. Con- 
sequently, SOC is only seen when a > 1. We see that 
the Liapunov exponents increase as a or Sf increases, 
with the other parameters kept fixed. Fig. ^(c) shows 
the largest Liapunov exponent as a function of the sys- 



tem size for a = 4 and Sf — 0.1 (solid), a — 4 and 
Sf = 0.01 (dashed) and a = 2.5 and Sf = 0.01 (short- 
dashed). We observe that Ai is approximately constant 
for small L and decreases nearly as power-law when L is 
greater than a certain value. The value where this bend- 
ing occurs seems to be sensitive to both Sf and a. Since 
Ai — > in the thermodynamic limit (L — > 00) we con- 
clude that chaos exists only in a microscopic level, and 
that any time or space scales in the system are negligi- 
ble. We have studied the system using slower time scales, 
such as measuring time by the updates in the first block. 
Still we find that Ai — > when L — > 00. In the train 
model governed by ODE's and pulled with a constant 
finite velocity we have found that the largest Liapunov 
exponent tends to a constant as the system size increases 
|l5| . However, our new unpublished results shows that 
the Liapunov exponents in that system start to decrease 
for system sizes greater than a given value, as it happens 
in the system we introduce here. 
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FIG. 1. (a) The train model, (b) The function <f>{x). 



FIG. 2. Probability distribution P(s) of the number of 
toppling events s. (a) solid: L = 512, a = 4, Sf = 0.1; 
dashed: L = 512, a = 4, Sf = 0.01; short-dashed: L = 512, 
a = 2.5, Sf = 0.01. (b) L = 32, 64, 128, 256, 512 with a = 4 
and Sf = 0.1. (c) Data collapse of the cases shown in (b) with 
r = 1.54 and D = 2.20. 
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FIG. 3. (a) Probability distribution P(T) of the number 
of toppling events T for L = 32, 64, 128, 256, 512 with a = 4 
and 5/ = 0.1. (b) Data collapse of the cases shown in (a) 
using P(T) = T-yfriT/L") with y = 1.84 and a = 1.40. 

FIG. 4. (a) The largest (solid) and the second largest 
(dashed) Liapunov exponents as a function of a for Sf = 0.1, 
and (b) as a function of Sf for a = 4, with L = 64 in both 
cases, (c) The largest Liapunov exponent as a function of L 
for a = 4 and 5/ = 0.01 (dashed), a = 2.5 and 5/ = 0.01 
(short-dashed). The error bars associated with this calcula- 
tion have approximately the size of the plus sign shown on 
the curves. 
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Fig. 1 
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Fig. 2(a) 
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Fig. 4(c) 
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